// Based on sum/gsl_sum.h
// Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
// Original Author:  G. Jungman
// Modified By: M. Jarvis


#ifndef SERIESBASE_H
#define SERIESBASE_H

#include <vector>
#include <complex>
#include <algorithm>
#include <ostream>

#include <assert.h>

#define INF 1.e500

// The specific cases: (Each in their own .h file)
// SeriesLevin uses Levin u transform 
// SeriesEpsilon uses Wynn's Epsilon transform
// SeriesTheta uses Brezinski's Theta transform
// SeriesAitken uses Aitken's iterated Delta^2 transform

template <class T>
class SeriesBase
{
#define INITALL \
    doerrcalc(_doerrcalc), saveterms(_saveterms), nterms(0), sum_plain(0.),\
    sum_var(0.), actual_trunc_n(0.0), trunc_n(0.0), noise_n(0.0),\
    converging(false), better_n(false), pastbest(false), \
    best_result(0.0), best_trunc(INF), best_noise(INF), best_n(0),\
    zero_terms(0)

#define INITVECS(n) \
    do { \
        if(saveterms) terms.reserve(n); \
        if (doerrcalc) termvars.reserve(n); \
    } while (false)

public:
    // Constructing Series with doerrcalc=false will execute a
    // quicker version which doesn't calculated the noise.
    // This is ok if you know how many terms will be needed (eg. if 
    // a similar sum has already been performed).
    // However, if array has too many terms, this version will 
    // actually be slower than the non-"Quick" version.
    SeriesBase(uint n, bool _doerrcalc, bool _saveterms) : INITALL
    { INITVECS(n); }
    virtual ~SeriesBase() {}

#undef INITALL
#undef INITVECS

    template <class U>
    void AddTerm (const T& term, const U& termerr);

    void AddTerm (const T& term) {AddTerm(term,fabs(term)*Epsilon());}

    void AddTerms (const std::vector<T>& v)
    { for(uint i=0;i<v.size();i++) AddTerm(v[i]); }

    template <class U>
    void AddTerms (const std::vector<T>& v, const std::vector<U>& verr)
    {
        assert(v.size() == verr.size());
        for(uint i=0;i<v.size();i++) AddTerm(v[i],verr[i]); 
    }

        template <class U>
        void Extrapolate(
            T* sum_accel, U* abserr, int* terms_used=0, std::ostream* fout=0);

        void Extrapolate (T* sum_accel)
        { Extrapolate(sum_accel,(double *)0); }

        uint NTerms() const 
        {
            if(saveterms) assert(nterms==terms.size());
            return nterms; 
        }
        const std::vector<T>& GetTerms() const 
        {
            assert(saveterms); return terms;
        }
        const T& PlainSum(double* plainsumerr=0) const 
        {
            if (plainsumerr) *plainsumerr = sqrt(sum_var);
            return sum_plain;
        }
        bool PastBest() const 
        { return pastbest; }

protected:
        // Actually calculate result and dsum
        virtual void CalcResult(const T& term)=0; 
        
        virtual const T& GetResult() const=0;
        virtual const std::vector<T>& GetDSum() const=0;

        // whether to check noise of calculation
        // (requires more memory, but stops calc sooner w/ better error estimate)
        const bool doerrcalc; 

        // whether to save the terms along the way
        const bool saveterms; 

private:
        std::vector<T> terms; // series terms
        std::vector<double> termvars; // variances of series terms
        uint nterms;

        double Noise() const 
        {
            if (!doerrcalc) return 0.;
            double variance = 0.;
            const std::vector<T>& dsum = GetDSum();
            for(uint i=0;i<termvars.size();i++) {
                variance += Norm(dsum[i])*termvars[i];
            }
            return sqrt(variance);
        }

        T sum_plain; // simple sum of series
        double sum_var;

        double actual_trunc_n;  // Truncation error for this result
        double trunc_n;  // Average of last two truncation errors
        double noise_n;  // Noise for this result

        bool converging; // Are we in the converging region?
        bool better_n; // Is this result better than the last?
        bool pastbest; // Are we past the good truncation region

        T best_result; // Extrapolation with smallest truncation error
        double best_trunc; // Truncation error for best_result
        double best_noise; // Noise for best_result
        int best_n; // N terms used for best_result

        int zero_terms;

        static double Epsilon();
        static double Norm(const T& x);
};


#define DBLEPS 2.2204460492503131e-16

template <> 
inline double SeriesBase<double>::Epsilon() {return DBLEPS;}

template <> 
inline double SeriesBase<std::complex<double> >::Epsilon() {return DBLEPS;}

template <> 
inline double SeriesBase<double>::Norm(const double& x) { return x*x; }

template <> 
inline double SeriesBase<std::complex<double> >::Norm(
    const std::complex<double>& x) 
{ return std::norm(x); }

#ifdef EXTRAPREC_H
#define QUADEPS 3.08148791094e-33
template <> 
inline double SeriesBase<Quad>::Epsilon() {return DBLEPS;}

template <> 
inline double SeriesBase<std::complex<Quad> >::Epsilon() {return DBLEPS;}

template <> 
inline double SeriesBase<Quad>::Norm(const Quad& x) { return x*x; }

template <> 
inline double SeriesBase<std::complex<Quad> >::Norm(const std::complex<Quad>& x) 
{ return std::norm(x); }
#endif

#if 0
template<class T> inline double SeriesBase<T>::Epsilon() 
{
    static double eps = 0.;
    if (eps == 0.0) {
        T x(1.);
        T teps(1.);
        while (true) {
            T y = x;
            y += teps;
            y -= x;
            if (fabs(y) < fabs(teps)/2.) {
                eps = double(real(teps));
                //std::cerr<<"epslon = "<<eps<<std::endl;
                break;
            }
            teps /= 2.;
            x += teps*teps;
        }
    }
    return eps;
}
#endif

template <class T> template <class U>
inline void SeriesBase<T>::Extrapolate (
    T* sum_accel, U* abserr, int* terms_used, std::ostream* fout)
{
    if(fout) (*fout)<<"In Extrapolate: \n";
    uint nterms = NTerms();
    if (nterms == 0) {
        if(fout) (*fout)<<
            "No (non-zero) terms in Series, "<<zero_terms<<" zero terms\n";
        if (zero_terms >= 2) {
            *sum_accel = T(0.); 
            if(abserr) *abserr = 0.0;
            if(terms_used) *terms_used = zero_terms;
        } else {
            *sum_accel = T(0.); 
            if(abserr) *abserr = INF;
            if(terms_used) *terms_used = zero_terms;
        }
    } else if (nterms == 1) {
        if(fout) (*fout)<<
            "Only 1 (non-zero) term in Series, "<<zero_terms<<" zero terms\n";
        if (zero_terms >= 2) {
            *sum_accel = terms[0]; 
            if(abserr) *abserr = 0.0;
            if(terms_used) *terms_used = 1 + zero_terms;
        } else {
            *sum_accel = terms[0]; 
            if(abserr) *abserr = INF;
            if(terms_used) *terms_used = 1 + zero_terms;
        }
    } else {
        assert(nterms > 1);
        if (fout) for(uint i=0;i<nterms;i++) {
            (*fout)<<i<<"  "<<terms[i]<<"  ";
            if (doerrcalc) {(*fout)<<sqrt(termvars[i]);}
            (*fout)<<std::endl;
        }

        if (converging) {
            // Stopped in the convergence region.  Return result and
            // error estimate.  
            *sum_accel = best_result;
            if(abserr) *abserr = std::max(Epsilon(),best_trunc + best_noise);
            if(terms_used) *terms_used = best_n+1;
            if(fout) (*fout)<<"Best = "<<best_result<<
                " +- "<<best_trunc+best_noise<<" using "<<best_n+1<<" terms\n";
        } else {
            // Never reached the convergence region.  Use the last
            // calculated values.  
            *sum_accel = GetResult();
            if(abserr) *abserr = trunc_n + noise_n;
            if(terms_used) *terms_used = nterms;
            if(fout) (*fout)<<"Not converging yet.... Sum = "<<
                GetResult()<<std::endl;
        }
    }
}

template <class T> template <class U>
inline void SeriesBase<T>::AddTerm(const T& term, const U& termerr)
{
    double termvar;
    if (doerrcalc) {
        double temp = std::max(double(termerr),Epsilon()*fabs(term));
        termvar = temp*temp;
    } else {
        termvar = 0.;
    }

    assert(!doerrcalc || !saveterms || termvars.size() == terms.size());

    int n=NTerms();

    // Ignore pass if already past best or this term is exactly 0
    if (pastbest) return;
    if (term==T(0.0)) { zero_terms++; return;}

    T result_nm1 = (n==0 ? T(0.) : GetResult());

    sum_plain += term;
    nterms++;
    if (saveterms) terms.push_back(term);
    if (doerrcalc) {
        sum_var += termvar;
        termvars.push_back(termvar);
    }
    CalcResult(term);

    T result_n = GetResult();

    // Check convergence properties
    const double SMALL = 0.01;

    // Compute the truncation error directly 
    double actual_trunc_nm1 = actual_trunc_n;
    actual_trunc_n = fabs (result_n - result_nm1);

    // Average results to make a more reliable estimate of the
    // real truncation error 
    double trunc_nm1 = trunc_n;
    trunc_n = 0.5 * (actual_trunc_n + actual_trunc_nm1);

    // Determine if we are in the convergence region.  
    bool better_nm1 = better_n;
    double absres = fabs(result_n);
    better_n = (trunc_n < trunc_nm1 || trunc_n < SMALL * absres);
    converging = converging || (better_n && better_nm1);

    if (converging) {
        //noise_n = std::min(Noise(),sqrt(sum_var));
        noise_n = Noise();
        if (trunc_n + noise_n < best_trunc + best_noise ) {
            // Found a low truncation point in the convergence region. Save it. 
            best_trunc = trunc_n;
            best_result = result_n;
            best_noise = noise_n;
            best_n = n;
        }
        if (noise_n > trunc_n && best_n != n) {
            pastbest = true;
        }
        if (trunc_n < 10.0 * Epsilon() * absres) {
            pastbest = true;
        }
    }
}

#undef INF

#endif
